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Abstract 

We  give  a  brief  review  of  hysteresis  in  viscoelastic  polymers.  The  efforts  surveyed 
range  from  phenomenological  to  molecular  modelling  with  applications  involving  recent 
efforts  on  elastomers  to  biotissue. 


1  Introduction 

In  control  and  systems  theory,  delay  or  hereditary  systems  or  systems  with  memory  (hys¬ 
teresis)  have  played  an  important  role  for  many  years  because  of  the  early  realizations  by 
Minorsky  and  others  [40,  41,  56,  57,  58,  60]  that  feedback  design  based  on  dynamics  wherein 
one  ignores  any  delays  may  fail  catastrophically  to  stabilize  or  control  a  system  in  which 
delays  or  hysteresis  are  present  in  the  dynamics.  This  is  true  whether  the  hysteresis  is  a 
fundamental  part  of  the  underlying  dynamics  or  a  part  of  the  input  or  control  operator. 
For  the  latter  there  is  a  growing  body  of  literature  [11,  12,  13,  45,  54,  68]  on  the  Preisach 
and  related  theories  for  hysteretic  control  input  such  as  arises  in  smart  material  systems 
[27,  39,  64],  Applications  in  which  delays  and  hysteresis  play  a  basic  role  in  the  underly¬ 
ing  dynamics  include  sustained  efforts  in  biology  with  early  efforts  employing  delay  systems 
[1,  2,  16,  17,  28,  29,  30,  31,  32,  42,  51,  53,  55,  62]  and  more  recent  investigations  involving 
hysteretic  probabilistic  structures  [5,  7]  as  well  as  classical  materials  and  electromagnetics 
research  (see  [8]  and  the  references  therein).  These  applications  drove  a  substantial  amount 
of  mathematical  and  computational  research  on  hysteretic  systems  in  the  last  half  of  the 
20th  century,  e.g.,  see  [26,  40,  41,  46,  47]  among  the  many  books,  research  monographs,  and 
research  articles  written.  Here  we  shall  focus  on  the  delays  or  hysteresis  arising  in  the  fun¬ 
damental  dynamics  of  the  systems  to  be  stabilized  or  controlled.  In  particular  we  consider 
viscoelastic  materials  that  are  polymeric  in  nature.  This  includes  a  wide  range  of  materials 
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of  current  importance  such  as  (rubber  or  silicone  based)  filled  elastomers  and  all  types  of 
biotissue  (soft  tissue,  ligaments,  cartilage,  etc.). 

2  Viscoelasticity 

The  mathematical  modelling  of  viscoelasticity  (sometimes  also  loosely  referred  to  as  hystere¬ 
sis)  in  materials  using  ideas  from  elasticity  has  attracted  the  attention  of  a  large  number  of 
investigators  over  the  past  century.  Among  significant  contributors  (see  the  many  references 
in  [33,  34,  35,  37,  38,  52,  59,  61,  65,  67,  69,  70])  have  been  some  of  the  true  giants  from  the 
fields  of  engineering  and  material  sciences.  One  of  the  most  widely  used  empirical  models  for 
viscoelasticity  in  materials  is  the  Boltzmann  convolution  law  [22,  35,  37,  38,  69],  one  form 
of  which  is  given  in  equation  (1) 

A  d 

a(t)  =  ge(e(t))  +  CDe(t)  +  /  Y{t  —  s)—gv  (e(s),  e(s))  ds,  (1) 

J- oo  ds 

where  e  is  the  infinitesimal  strain,  Y  is  the  convolution  memory  kernel,  and  ge  and  gv  are 
nonlinear  functions  accounting  for  the  elastic  and  viscoelastic  responses  of  the  elastomers, 
respectively;  for  summaries  and  further  references,  see  Chapter  2  of  [38]  as  well  as  [22] .  This 
form  of  model,  when  incorporated  into  force  balance  laws,  results  in  integro-partial  differential 
equations  which  are  most  often  phenomenological  in  nature  as  well  as  being  computationally 
challenging  both  in  simulation  and  control  design.  This  stress-strain  law  implies  that  the 
stress  depends  not  only  on  the  current  strain  and  strain  rate  but  also  on  the  history  of 
the  strain  and  the  strain-rate.  It  is  very  important  to  note  that  the  stress-strain  law  (1) 
contains  various  standard  internal  strain  or  internal  variable  formulations  as  special  cases. 
The  anclastic  displacement  field  (ADF)  models  of  Lesientre  [48,  49]  for  composite  materials 
exhibiting  both  elastic  and  anelastic  displacement  fields  are  formulated  on  the  assumption 
that  the  host  elastic  material  contains  anelastic  materials  with  internal  strains  £\  which  are 
clastic  strain  driven.  That  is,  the  constitutive  laws  have  the  form 

a(t)  =  Ee(t )  -  Ai£i(t),  (2) 


where  the  internal  strain  is  given  by 

^i(^)  "1 — £i(^)  =  c2 £(t)i  £i(0)  —  0,  (3) 

r 

or  equivalently, 


£\{t)  =  /  c2e  T  s(s)ds. 

Jo 

Several  generalizations  of  this  formulation  exist,  e.g.,  Johnson,  et  al.,  [43,  44],  suggest  that 
the  internal  strain  is  strain  rate  driven,  i.e., 


£i  (t)  +  -£i(i)  =  c2e(t). 

T 


(4) 
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The  Boltzmann-type  law  (1)  (under  appropriate  assumptions  on  the  past  memory  from  —  oo 
to  0)  corresponds  to  an  internal  strain  model  of  the  form 

£i(t)  +  -ei(t)  =  ^i9v{e{t),e{t)),  ei(0)  =  0.  (5) 

T  QX 

This  form  is  often  chosen  since  one  finds  that  neither  (3)  nor  (4)  provide  laws  that  readily 
describe  experimental  data,  especially  in  the  cases  of  filled  elastomers,  biotissues  and  other 
molecular  polymers. 

3  Fung’s  Quasi-linear  Models 

Fung,  in  his  extensive  efforts  [38]  with  biomechanics  and  biotissue,  develops  and  presents  the 
quasi-linear  viscoelastic  constitutive  equation 

Sij(t)  =  Jf  Gijkiit  -  r)dS<kl)^T)]dr,  (6) 

where  StJ  is  the  Kirchoff  stress  tensor,  E  is  the  Green’s  strain  tensor,  Gijki  is  a  reduced 
relaxation  function,  and  Sff  is  the  “elastic”  stress  tensor.  For  the  scalar  components  Gijki, 
Fung  proposes  the  reduced  relaxation  function  G(t)  given  in  the  form 

G(t)  =  {l  +  C[Ei{— )  -  Ei(  —  )]|  [1  +  cln(  — )]  1.  (7) 

L  T2  T\  J  T\ 

Here  Ei(z)  =  f°°  C  represents  the  degree  to  which  viscous  effects  are  present,  and  T\ 

and  T2  represent  fast  and  slow  viscous  time  phenomena.  We  note  that  the  internal  strain 
variable  formulation  (2),  (5)  is  equivalent  to  the  constitutive  relationship  proposed  by  Fung 
if  one  considers  an  approximation  of  the  relaxation  function  G  by  a  sum  of  exponential 
terms.  Various  internal  strain  variable  models  are  investigated  in  [3]  and  a  good  agreement 
is  demonstrated  between  a  two  internal  strain  variable  model  (e.g.,  of  the  form  a  =  Ee  — 
EiEi  —  E2£2)  and  undamped  simulated  data  based  on  the  Fung  kernel  G. 

Since  its  introduction,  this  quasi-linear  viscoelastic  (QLV)  theory  of  Fung  has  been  ap¬ 
plied  successfully  in  stress-strain  experiments  to  several  types  of  biological  tissue.  A  benefit 
to  using  (6)  as  a  constitutive  equation  is  that,  unlike  simpler  models  for  viscoelasticity,  it 
allows  for  the  consideration  of  a  continuous  spectrum  (e.g.,  see  the  discussions  in  [38])  of 
relaxation  times  and  frequencies  (this  is  also  true  of  the  probabilistic-based  internal  variable 
approach  developed  in  [23]  and  described  below).  (The  need  for  a  continuum  of  relaxation 
times  in  certain  materials  was  observed  many  years  ago  [36,  63,  66,  70].)  While  Fung’s 
theory  has  been  successfully  employed  for  fitting  hysteretic  stress-strain  curves,  for  control 
applications  one  is  interested  in  using  it  in  a  full  dynamical  model.  Unfortunately,  the  QLV, 
as  presented  by  Fung,  leads  to  exceedingly  difficult  computations  within  full  dynamical  par¬ 
tial  differential  equations,  especially  in  estimation  and  control  problems.  This  motivated 
the  development  of  the  internal  variable  approach  described  in  [3,  23,  48]  (which  permits 
discrete  approximation  to  a  continuum)  in  attempts  to  approximate  well  the  correspond¬ 
ing  dynamic  responses  even  in  cases  where  the  stress-strain  curves  alone  do  not  produce 
adequate  approximations  -  see  [38]. 
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The  probabilistic  based  internal  variable  alternative  [23]  to  Fung’s  kernel  involves  a  pa¬ 
rameter  dependent  kernel  with  a  continuous  distribution  of  parameters  and  internal  variables. 
In  the  case  of  a  finite  combination  of  Dirac  6  distributions,  one  obtains  a  finite  summation  of 
exponential  functions  as  the  approximation  kernel  (see  the  discussions  below).  This  method 
can  be  extended  to  allow  for  consideration  of  a  continuous  spectrum  of  relaxation  times 
and  frequencies  by  utilizing  absolutely  continuous  parameter  distributions  in  place  of  the  5 
distributions. 


4  Internal  Variables 

The  internal  variable  approach  to  overcome  both  conceptual  and  computational  challenges  is 
consistent  with  the  belief  that  hysteresis  is  actually  a  manifestation  of  the  presence  of  multi¬ 
ple  scales  in  a  physical  or  biological  material  system  that  is  frequently  modelled  (and  masked) 
with  a  phenomenological  representation  such  as  an  hysteresis  integral  for  the  macroscopic 
stress-strain  constitutive  law.  The  internal  variable  modelling  leads  to  an  efficient  compu¬ 
tational  alternative  for  the  corresponding  integro-partial  differential  equation  models.  In 
addition,  it  provides  a  “molecular”  basis  for  the  models  (for  a  comparison  of  models  of  vis¬ 
coelastic  damping  via  hysteretic  integrals  versus  internal  variable  representations,  see  [22] 
and  the  references  therein). 

Our  own  interest  in  viscoelasticity  in  polymeric  materials  has  been  motivated  by  projects 
in  our  Industrial  Applied  Mathematics  Program  with  at  least  two  of  our  industrial  partners: 
The  Lord  Corporation  and  Medacoustics,  Inc.  The  collaborations  with  polymer  scientists 
and  engineers  at  Lord  involved  the  dynamic  modelling  of  filled  rubbers  which  experimen¬ 
tally  exhibit  both  significant  hysteresis  and  nonlinearity  in  tensile  and  shear  deformations  as 
depicted  in  the  sample  stress-strain  curves  in  Figure  1.  The  efforts  with  engineers  at  Meda¬ 
coustics  used  some  of  the  viscoelastic  models  we  have  investigated  in  attempts  to  understand 
the  propagation  of  arterial  stenosis  induced  shear  waves  in  composite  biotissue  in  a  sensor 
development  and  characterization  project. 


Figure  1:  Experimental  stress-strain  curves  for  (1)  unfilled,  (2)  lightly  filled  and  (3)  highly 
filled  rubber  in  tensile  deformations. 


4 


In  some  of  our  earlier  efforts  [24,  25],  the  models  for  hysteretic  damping  in  elastomers 
employed  a  phenomenological  Boltzmann-type  constitutive  law  of  the  form  (1).  As  explained 
in  [20,  24],  our  nonlinear  materials  undergoing  large  deformations  required  the  use  of  finite 
(as  opposed  to  infinitesimal)  strain  theories  [61].  However,  since  the  nonlinearity  between  the 
stress  and  finite  strain  is  an  unknown  to  be  estimated  (using  inverse  problem  algorithms)  and 
since  the  finite  strain  can  be  expressed  in  terms  of  known  nonlinearities  as  a  function  of  the 
infinitesimal  strain  (at  least  in  the  problems  of  interest  here),  one  can  effectively  formulate 
the  problem  as  one  of  estimating  the  unknown  nonlinearity  between  stress  and  infinitesimal 
strain  (see  [24]).  Hence  one  can  develop  models  for  stress  in  terms  of  infinitesimal  strain.  Our 
previous  efforts  as  summarized  in  [20,  21]  have  shown,  through  comparison  with  experimental 
data,  that  the  best  fit  to  filled  elastomer  data  occurs  when  ge  and  gv  are  cubic,  along  with  Y  as 
a  distribution  of  decaying  exponentials.  We  subsequently  [18,  19,  20,  21]  developed  nonlinear 
models  based  on  stick-slip  “molecular”  ideas  of  Johnson  and  Stacer  [43]  and  Doi  and  Edwards 
[34]  which  resulted  in  a  form  for  ge,  gv  and  Y  in  (1)  that  matched  the  empirical  hirelings 
reported  in  [20,  24,  25].  These  models  allow  for  multiple  relaxation  times  present  in  polymer 
strands  of  composite  materials  within  a  virtual  compartmental  model  of  entangled  chemically 
cross-linked/physically  constrained  system  of  long  chain  “molecules”.  While  accounting  for 
multiple  relaxation  parameters,  these  models  do  not  include  physically  or  chemically  based 
parameters  in  representations  of  the  polymer  strands. 

We  mention  briefly  two  recent  advances:  ( i )  a  new  constitutive  model  [10]  that  has 
been  developed  which  combines  the  virtual  stick-slip  continuum  “molecular-based”  ideas  of 
Johnson  and  Stacer  [43]  with  the  Rouse  bead  chain  (see  Figure  2  below)  ideas  as  described 
in  Doi  and  Edwards  [34];  {ii)  a  two  dimensional  version  [14,  15,  50]  of  a  model  that  accounts 
for  stenosis  driven  shear  wave  propagation  in  biotissue. 


5  A  Stick-Slip/Rouse  Hybrid  Model 

The  early  models  and  the  nonlinear  extensions  of  the  Boltzmann  law  did  not  provide  insight 
into  the  underlying  mechanisms  for  tensile  and/or  shear  deformations  in  filled  rubber  or 
biotissue.  This  is  not  unexpected  since  the  approaches  described  above  are  based  on  pseudo- 
phenomenological  formulations.  We  then  [19,  20,  21]  turned  to  a  different  approach  based 
on  molecular  arguments  which,  as  we  shall  see,  lead  precisely  to  the  class  of  models  based  on 
a  Boltzmann  hysteresis  formulation.  As  usual,  one  begins  with  force  and  moment  balance 
and  seeks  constitutive  laws  for  the  viscoelastic  stress  term  crvisco  in 

® elastic i  T  ^visco(^l  ('))  i 

where  £  —  is  the  infinitesimal  strain  and  £\  is  an  “internal  strain”  variable  on  which  oxnsco 
depends  in  an  hysteretic  manner.  As  described  above,  we  found  that  a  reasonable  description 
of  the  data  of  interest  could  be  given  with  the  typical  stress-strain  relationship 

rt  t  d 

<r(t)  =  ffe(e(t),e(t))  +  J  'ye^  —gv{e{s),e{s))ds,  (8) 
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where  r  is  a  relaxation  parameter,  gv  is  defined  with  cubic  polynomials  and  ge  =  ge ( )  + 
AVe  haye  already  observed  that  this  expression  is  equivalent  to 

a(t)  =  ge(e(t),e(t))  +  je^t;  r), 

where,  for  a  given  “relaxation  parameter”  r,  the  internal  strain  £i(t',r)  satisfies  (5).  In 
fact,  we  found  that  highly  filled  rubbers  required  multiple  relaxation  times  Ti,T2  in  an  ap¬ 
proximation  to  (7)  to  obtain  good  model  fits  to  the  data.  As  one  might  expect,  molecular 
based  formulations,  where  microscopic  relaxation  parameters  vary  across  the  population  of 
molecules  in  the  material,  lead  to  internal  dynamics  of  the  form  (5), (8)  that  involve  multiple 
values  of  r.  When  combined  with  a  Prohorov  metric  framework  (see  [4,  6])  for  uncertainty  in 
internal  dynamics,  these  ideas  lead  to  the  computational  models  we  have  used.  Indeed,  the 
molecular  based  approach  leads  to  a  general  class  of  models  with  uncertainty  or  randomness 
in  the  stress 


a(t,x;P)  =  ge(e(t,x),e(t,x))+'y  j  £i (f,  a;;  r)cLP(r),  (9) 

where  P  is  a  probability  distribution  over  the  set  T  of  possible  relaxation  parameters,  and 
£i(t;  t )  satisfies,  for  each  r  G  T , 

£\ (t,  x ;  r)  H — £\ (t,  x ;  r)  =  i(t,  x)h(e(t,  x)). 

T 

For  the  reptation  model  derivation  in  [19],  one  begins  with  the  Doi/Edwards  [34]  stick-slip 
molecular  models  as  embodied  in  the  continuous  tube  reptation  models  of  Johnson/Stacer 
[43]  wherein  polymer  materials  such  as  rubber  are  postulated  to  be  composed  of  two  types  of 
molecules.  In  tensile  deformations,  one  denotes  by  L(t)  the  length  of  chemically  cross-linked 
or  CC  molecules,  while  £(t)  denotes  the  length  of  physically  constrained  or  PC  molecules. 
To  use  stick-slip  models  in  continuum  simulations  of  reptation  in  rubbers,  one  considers 
networks  of  “cells”  or  boxes  of  parallel-sided  CC  boxes  and  PC  boxes  with  sides  of  length 
(principal  stretches) 

,  .  .  duc  .  „  „  duv 

Ac  =  l  +  e  =  l  +  -7— ,  Ap  =  l  +  £i  =  l  +  — — , 

ox  ox 

respectively.  Here  uc  denotes  the  deformations  of  the  CC  box  and  up  denotes  the  defor¬ 
mations  of  the  PC  box.  Using  a  linear  stick-slip  assumption  as  in  [43],  and  strain  energy 
densities  based  on  experiments  of  Young  and  Danik  (see  [18,  20]  for  details),  one  obtains  as 
a  limit  of  PC  response  to  step  tensile  deformations  of  the  CC  molecules,  the  £,  £1  coupled 
dynamics 

.  ^  1  _  .1  +  ei 

r  1  +  £ 

However,  if  one  replaces  the  linear  assumption  of  [18]  by  a  nonlinear  stick-slip  hypothesis 
(which  is  the  basis  of  the  work  in  [19]),  one  obtains  a  more  general  nonlinear,  dynamical 
relationship  between  £  and  £1  given  by 


£l  H - £1  —  £/((  1  +  £i)/(l  +£))• 

T 
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Expansion  and  truncation  of  higher  order  terms  lead  to  equations  of  the  form 

X 

Sl  H - Si  =  s(cx o  +  OL\S  +  +  ci^s3),  (10) 

T 

which  are  of  the  same  form  as  the  internal  variable  model  (5), (8)  with  gv  a  cubic  polynomial. 
For  the  corresponding  contributions  to  a  from  the  strain  energy  densities  of  Young-Danik/ 
Johnson-Stacer  with  the  nonlinear  stick-slip  hypothesis,  one  obtains  a  contribution  to  the 
rate  independent  strain  g *  (after  expanding  /  in  a  Taylor  series  and  dropping  higher  order 
terms)  of  the  form 

9v(£i£  l)  —  9cubic{£ )  +  7l£l ) 

where  s i  is  as  before  (i.e.,  the  internal  strain  satisfying  (10)).  Thus,  the  total  stress-strain 
relationship  can  be  written  in  the  form  (9).  If  the  measure  P  of  (9)  has  atoms  at  tx  and 
r2,  (i.e.,  the  measure  is  composed  of  Dirac  measures  concentrated  at  tx  and  r2),  then  the 
constitutive  law  leads  precisely  to  the  model 

a(t,x]P)  =  ge(s(t,x),s(t,x ))  +  jiSi(t,  x;  n)  +  72£2(t,  x;  r2), 

which  was  used  in  the  data  fits  in  [3,  20].  Further  details  on  this  modeling  approach  can  be 
found  in  [9,  21] 

We  next  developed  a  new  constitutive  model  which  combines  the  “molecular-based” 
ideas  of  Johnson  and  Stacer  with  the  Rouse  bead  chain  ideas  and  explained  its  relation  to 
the  Boltzmann  phenomenological  models. 


Figure  2:  Representation  of  vectors  for  a  bead-spring  polymer  molecule. 

The  new  molecular-based  constitutive  model,  in  which  polymer  chains  are  treated  as 
Rouse  type  strings  of  interconnected  beads  (a  reasonable  approximation  for  many  materi¬ 
als),  permits  the  incorporation  of  many  important  physical  parameters  (such  as  temperature, 
segment  bond  length,  internal  friction,  and  segment  density)  in  the  overall  hysteretic  consti¬ 
tutive  relationship.  Its  form  is  similar  to  that  developed  in  [19,  20]  and  does  have  the  general 
form  (1)  of  Boltzmann  type,  even  though  the  kernel  is  not  of  convolution  type.  The  result¬ 
ing  model  does,  however,  provide  a  molecular  basis  for  the  earlier  pseudo-phenomenological 
stick-slip  models. 

We  give  only  a  brief  outline  of  the  new  constitutive  model  here;  more  details  of  the  deriva¬ 
tion  can  be  found  in  the  report  [10].  We  model  a  polymer  material  undergoing  directional 
deformation  by  assuming  it  is  composed  of  two  virtual  compartments  as  depicted  in  Figure  3. 
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Entrapped  Portion  of  PC  molecule 


Escaped  Portion  of  PC  molecule 


Figure  3:  PC  molecule  entrapped  by  the  surrounding  constraining  tube. 

One  compartment  consists  of  a  constraining  tube  which  is  a  macroscopic  compartment  con¬ 
taining  both  CC  (chemically  cross-linked)  and  PC  (physically  constrained)  molecules.  The 
other  compartment  is  microscopic  in  nature  and  consist  of  those  PC  molecules  aligned  with 
the  direction  of  the  deformation.  These  molecules  will  at  first  “stick”  to  the  constraining 
tube  and  be  carried  along  with  its  motion,  but  will  very  quickly  “slip”  and  begin  to  “relax” 
back  to  a  configuration  of  lower  strain  energy.  In  the  model  derivation  one  computes  the  con¬ 
tributions  of  both  “compartments”  to  the  overall  stress  of  this  polymer  material  undergoing 
deformations  to  obtain  the  constitutive  law. 


6  Stenosis-Driven  Shear  Wave  Propagation  in  Biotis¬ 
sue 

Finally,  we  turn  to  recent  results  on  the  viscoelastic  models  for  propagation  of  stenosis-driven 
biotissue  waves  mentioned  above.  In  the  biotissue  efforts  [14,  15,  50],  we  employ  an  internal 
variable  formulation  of  Boltzmann  type  hysteresis  laws  to  investigate  the  propagation  of 
stenosis  generated  waves  in  biotissue  where  it  has  been  demonstrated  that  a  viscoelastic 
(as  opposed  to  an  elastic)  formulation  is  important  and  that  waves  generated  in  a  two- 
dimensional  cylindrical  geometry  with  inner  radius  partial  occlusions  can  be  readily  modelled 
and  simulated. 

Specifically  we  have  carried  out  efforts  on  two  and  three  dimensional  models  that  employ 
an  internal  variable  approach  to  model  wave  propagation.  To  motivate  this,  we  recall  [3] 
that  coronary  artery  disease  (CAD)  is  caused  by  atherosclerosis,  the  gradual  accumulation 
of  plaque  along  the  walls  of  an  artery.  This  buildup,  known  as  a  stenosis,  restricts  the  flow  of 
blood,  leading  to  a  decrease  in  the  oxygen  supply  to  the  heart  muscle.  It  is  well  known  that 


arterial  stenoses  produce  sounds  due  to  turbulent  blood  flow  in  partially  occluded  arteries. 
In  principle,  turbulent  normal  wall  forces  exist  at  and  downstream  from  an  arterial  stenosis, 
exerting  pressure  on  the  wall  of  the  artery  which  then  causes  a  small  displacement  in  the 
surrounding  body  tissue.  The  goal  is  to  model  the  propagation  of  the  wave  generated  from 
the  stenosis  to  the  chest  wall,  and  ultimately,  to  create  an  inverse  problem  methodology 
which  can  be  utilized  to  determine  the  location  of  an  arterial  stenosis.  In  [14,  15,  50]  we 
also  discuss  comparison  of  the  viscoelastic  model  to  an  clastic  one  as  well  as  present  typical 
simulations  for  a  biologically  motivated  example. 
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